Disorder Induced Localized States in Graphene 
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We consider the electronic structure near vacancies in the half-filled honeycomb lattice. It is shown that 
vacancies induce the formation of localized states. When particle-hole symmetry is broken, localized states 
become resonances close to the Fermi level. We also study the problem of a finite density of vacancies, obtaining 
the electronic density of states, and discussing the issue of electronic localization in these systems. Our results 
have also relevance for the problem of disorder in d-wave superconductors. 
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Introduction. The problem of disorder in systems with 
Dirac fermions has been studied extensively in the last few 
years in the context of dirty d-wave superconductors [1]. 
Dirac fermions are also the elementary excitations of the hon- 
eycomb lattice at half-filling, equally known as graphene, 
which is realized in two-dimensional (2D) Carbon based ma- 
terials with sp 2 bonding. It is well-known that disorder is 
ubiquitous in graphene and graphite (which is produced by 
stacking graphene sheets) and its effect on the electronic struc- 
ture has been studied extensively [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
12, 13]. It has been shown recently [14] that the interplay 
of disorder and electron-electron interactions is fundamental 
for the understanding of recent experiments in graphene de- 
vices [15]. Furthermore, experiments reveal that ferromag- 
netism is generated in heavily disordered graphite samples 
[16, 17, 18, 19, 20], but the understanding of the interplay 
of strong disorder and electron-electron interactions in these 
systems is still in its infancy. Different mechanisms for ferro- 
magnetism in graphite have been proposed and they are either 
based on the nucleation of ferromagnetism around extended 
defects such as edges [2, 8, 12, 13] or due to exchange in- 
teractions originating from unscreened Coulomb interactions 
[21]. Therefore, the understanding of the nature of the elec- 
tronic states in Dirac fermion systems with strong disorder is 
of the utmost interest. 

In the following, we analyze in detail states near the Fermi 
energy induced by vacancies in a tight-binding model for the 
electronic states of graphene planes. We show that single va- 
cancies in a graphene plane generate localized states which 
are sensitive to the presence of particle-hole symmetry break- 
ing. Moreover, a finite density of such defects leads to strong 
changes in the local and averaged electronic Density Of States 
(DOS) with the creation of localized states at the Dirac point. 

The model. We consider a single band model described by 
the Hamiltonian: 
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where Cj (cj) annihilates (creates) an electron at site Rj of 
the honeycomb lattice (the spin quantum numbers are sup- 
pressed since we do not consider spin dependent phenomena). 



In (1) t is the nearest neighbor hopping energy (t sa 2.7eV in 
graphene) and tf is the next-nearest neighbor hopping energy. 
In the honeycomb lattice, while t describes the hopping of 
electrons between the two sub-lattices, if describes the hop- 
ping in the same sub-lattice and therefore breaks particle-hole 
symmetry. The relevance of t' for graphene is suggested by 
different experiments (t' w 0.2f) [14]. 

Localized states around vacancies. Firstly we consider the 
particle-hole symmetric case (t 1 = 0), and use the geometry 
shown in Fig. 1 . We analyze a cluster with periodic bound- 
ary conditions along the vertical direction, like a zigzag nan- 
otube. In the absence of vacancies, the Hamiltonian can be 
simplified using the translational symmetry along the verti- 
cal direction. The states can be classified by the momentum 
(in units of 1/a, with a the lattice constant) along the ver- 
tical axis, k m = {2itm)/N , m = 1,...,N where N is the 
number of unit cells along the vertical axis. The system is 
metallic if N = 3M, where M is an integer. Using this 
set of momenta, the wavefunction amplitudes in each sub- 
lattice (A,B) can be written as aij — J2k a z.fc m e ifcm "' and 
— J2k ^iM m & %km K where are the unit cell coor- 
dinates in the geometry of Fig. 1. Upon this transformation, 
the original problem is mapped, for each value of k m , into a 
one-dimensional (ID) tight-binding model along the horizon- 
tal direction, with two sites per unit cell. Such ID problem is 
characterized by site amplitudes a/ = ai t k m and 6/ = bi ; k m , 
and by two effective hoppings, t and t(l + e lkm ); a gauge 
transformation allows us to make them both real, reading t 
and 2t cos(fc m /2). 

For the solution of the impurity problem, it is convenient 
to define two planes contiguous to the vacancy, as shown in 
Fig. 1 . The two planes belong to the same sub-lattice, opposite 
to the one where the vacancy resides. A possible localized 
state at zero energy must: (/) decay as one moves along the 
horizontal axis away from these planes, (ii) satisfy the bulk 
tight binding equations arising from the eigenvalue condition 
H\ty) = E\*i?) = at these axes and beyond, and (Hi) have 



amplitudes a^' and a^' which satisfy the equations, 
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FIG. 1 : (color online) Sketch of the geometry considered in the text 
for the study of a single B-site vacancy (®). The system has periodic 
boundary conditions along the vertical axis, and is infinite along the 
horizontal axis. Only some A-site amplitudes are shown. 



everywhere automatically satisfies condition (;;). Plane L at 
the left of the vacancy can be considered a zig-zag edge. To 
its left, one can define states which decay exponentially for 
momenta such that 27r/3 < k m < Att/3 [5, 6]. The asso- 
ciated wavefunctions decay as |2cos(fc rn /2)|'. As defined in 
Fig. 1, / is always positive, since it grows away from both the 
L and the R planes, and can be interpreted as the distance to 
plane L divided by 3a/2. Analogously, one can define local- 
ized wavefunctions to the right of plane R, for momenta such 
that < k m < 27r/3 or 47r/3 < k m < 2tt. These wave- 
functions decay as |2 cos(/s m /2)|~'. The amplitudes and 

ciQ R j at planes L and R can be written, in terms of momentum 
eigenstates, as: 

k m k m , 

and if there was no impurity present in the system, we would 
have k m and fe m ' m the range [0, 27r[. The boundary condition 
introduced by Eq. (2) can be written as: 

e ^ ikmj = -e (i + eife '-') a iy km,j . w 

for j 7^ 0. These equations admit the solution: 

a[ L 2 = l, a£(l + e^)=l» (5) 

where the wavevectors k m and k m > satisfy 27r/3 < k m < 
4?r/3 and < k m > < 2tt/3 or 4tt/3 < k m , < 2tt. We 
note that these momentum values are the same defining the 



decaying states to the left and to the right of planes L and R, 
respectively, if we had considered the two cases as separate 
problems. Hence, we can use the amplitudes we found in the 
latter case (given above) when constructing the wave function 
for an impurity. To the left of plane L, the amplitude a^ L ) is 

now given as a)*/") = J2k m [~ 2 cos(fc/2)f exp[ik m (j + 1/2)}, 
where the values of k m are those imposed by the boundary 
condition, and a similar expression for af^ . For sufficiently 
wide nanotubes (corresponding to the solution for the infi- 
nite lattice) we can approximate the sums in k m by integrals. 
Shifting from lattice position coordinates to distances rela- 
tive to the L plane, the amplitude gives the wavefunc- 

tion \I>(x = I3a/2,y = a\/3(j + 1/2)) at a point of coordi- 
nates (x, y). In units of the lattice constant, \&(x, y) is approx- 
imately 

/.4tt/3 

^ {L) (x,y) - / dk(-2cos(k/2)) 2x/3 e lky/ ^ 

J2it/3 

(47riy)/(3\/3) e 2wi(x+y/ V3)/3 

« —. + : , (6) 

x + ly x — ly 

when the lattice site (x, y) is in the opposite sub-lattice of the 
vacancy, and ^(x, y) = when (x,y) is in the same sub- 
lattice as the vacancy. It is now clear that extra vacancies 
added to the sub-lattice where the impurity resides, cannot 
change this wave function. 

We would like to understand this results, from the point of 
view of a low-energy effective theory. The eigenstates of the 
discrete Hamiltonian (l) can be approximated, at long wave- 
lengths, by the Dirac equation. The wavefunctions can be 
written as a spinor: 

where the functions ipi(x, y) and ip2(x, y) correspond to the 
amplitudes of the wavefunctions in each of the two sub- 
lattices. There are two set of spinors, a, b, which correspond 
to the two inequivalent states at the corners of the Brillouin 
zone. At zero energy, the functions at one of the Dirac points 
satisfy: d z tfji(z, z) = , dzip%(z, z) = 0, where z = x + iy 
and z = x — iy; and those at the other Dirac point can be ob- 
tained by replacing zby z everywhere. The result of Eq. (6) 
implies that the boundary conditions at the vacancy are such 
that the combination 

•(.,»)- (*"J* tw )«(ij*) w 

is selected (notice that, in the long wavelength limit, the phase 
factors in Eq. (6) are factored out upon defining the slowly 
varying Dirac fields). This solution, although decaying away 
from the impurity, is not normalizable. It has the same spa- 
tial dependence as the quasi-localized solutions which are in- 
duced by radial potentials on 2D Dirac fermions [22]. The 
matching of localized states described above cannot be gener- 
alized to the case t 1 ^ 0, as the band of edge states is not de- 
generate in energy [23]. The localized state at the Fermi level 
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FIG. 2: (color online) (a) Comparison between the local DOS 
in the vicinity of a vacancy (blue/continuous) with the bulk DOS 
(red/dashed) in clean systems, (b) Total DOS in the vicinity of the 
Dirac points for clusters with 4 x 10 6 sites, at selected vacancy con- 
centrations. Numerical results in (a,b) obtained for t' = (top pan- 
els), t' = 0.1 1 (center) and t' = 0.2 1 (bottom). (Notice the scale 
truncation in the upper part of the first panel in (a). 

becomes a resonance inside the continuum of extended states. 
Numerical results for the local (DOS) at a site near a single 
vacancy are shown in Fig. 2(a). In the absence of electron- 
hole symmetry, the localized state becomes a resonance with 
increasing width, shifted from the Fermi energy. 

The scheme used here can be generalized to study other 
lattice representations of the Dirac equation with electron- 
hole symmetry. The electronic structure near impurities in 
d-wave superconductors has been studied using the Hamilto- 
nian: [24, 25, 26]: 

Hbcs = t £ c\ Cj +A (cj jC ] ±aJ - ct . c ] . ±b ) + h.c. , (9) 

where the sites are defined in a square lattice, a and b are 
the lattice vectors along the horizontal and vertical directions. 
The Hamiltonian (9) is formally identical to a tight-binding 
model on a square lattice and two orbitals per site. The 
hopping terms between different orbitals have different signs 
along the horizontal and vertical axis. If periodic boundary 
conditions are applied along the (1,1) and (1,-1) directions, 
the problem can be written as the matching of a set of ID 
wavefunctions, in a way similar to the scheme depicted in 
Fig. 1. A vacancy leads to a quasi-localized state, described 
by the same long wavelength wavefunction (see also [27]). 
The quasi-localized state described here, being localized in a 
single sub-lattice, remains a solution when there are two va- 
cancies in different sub-lattices, in agreement with numerical 
calculations [28]. 

Finite densities of vacancies. We have extended the previ- 
ous results to systems with a finite density of vacancies, using 
the stochastic recursion method to obtain the DOS in clusters 



with up to ~ 10 6 sites. Results for different impurity con- 
centrations and different values of t' are shown in Fig. 2(b). 
In the presence of electron-hole symmetry (t' = 0), the in- 
clusion of vacancies brings an increase of spectral weight to 
the surroundings of the Dirac point, leading to a DOS whose 
behavior for E w mostly resembles the results obtained 
within a Coherent Potential Approximation (CPA) [14]. The 
most important feature, however, is the emergence of a sharp 
peak at the Fermi level, superimposed upon the fiat portion 
of the DOS (apart from the peak, the DOS flattens out in this 
neighborhood as x is increased past the 5 % shown here). The 
breaking of the particle-hole symmetry by a finite t' results 
in the broadening of the peak at the Fermi energy, and the 
displacement of its position by an amount of the order of t'. 
All these effects take place close to the the Fermi energy. At 
higher energies, the only deviations from the DOS of a clean 
system are the softening of the van Hove singularities and the 
development of Lifshitz tails (not-shown) at the band edge, 
both induced by the increasing disorder caused by the random 
dilution. The onset of this high energy regime, where the pro- 
file of the DOS is essentially unperturbed by the presence of 

— 1/2 

vacancies, is determined by e m Up /I, I ~ n i m p being the 
average distance between impurities. 

To address the degree of localization for the states near the 
Fermi level, the Inverse Participation Ratio (IPR) was cal- 
culated via exact diagonalization on smaller systems. For 
an eigenstate to, the IPR is the quantity defined as: V m — 
|^m(i) | 4 > the index i labeling the lattice sites. The wave- 
function of an extended state has an amplitude equally signifi- 
cant throughout the entire system (\l/ m (z) ~ N^ 1 ^ 2 ), whence 
we naturally expect V m ~ AT -1 . For a localized state, in 
opposition, its very definition entails the fact that only a finite 
number of lattice sites will contribute to the normalization, re- 
sulting in much higher values of V m ■ Results for different val- 
ues of t' are shown in Fig. 3 for random dilution at 0.5 %, One 
observes, first, that V m ~ 3 /AT for all energies but the Fermi 
level neighborhood, as expected for states extended up to the 
length scale of the system sizes used in the numerics. Sec- 
ondly, the IPR becomes significant exactly in the same energy 
range where the DOS exhibits the vacancy-induced anomalies 
discussed above. Clearly, the farther the system is driven from 
the particle-hole symmetric case, the weaker the localization 
effect, as illustrated by the results obtained with t' = 0.2 1. 
To this respect, it is worth mentioning that the magnitude of 
the strongest peaks in V m at t! =0 and t' = 0.1 1 is equal 
to the magnitude of the IPR calculated for a single impurity 
problem [29]. Such results indicate the existence of quasi- 
localized states at the center of the resonance, induced by the 
presence of the vacancies. For higher doping strengths, the 
enhancement of V m is weaker in the regions where the DOS 
becomes flat, in agreement with the reasonably good descrip- 
tion obtained with the CPA. 

Conclusions. We have studied the local DOS near vacan- 
cies in graphene planes. The global and local DOS, are valu- 
able tools in the connection between a microscopic theory and 
the interpretation of local spectroscopic experiments, of which 
recent scanning tunneling spectroscopy measurements are an 
example [20] . In agreement with general arguments valid for 
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finite concentration of impurities, we find a sharp peak built 
up from localized states, superimposed on a structureless fi- 
nite DOS, also induced by the vacancies. We have not consid- 
ered here additional potentials or local atomic displacements 
which can take place near vacancies. Our findings indicate 
that the main features in the electronic structure are extended 
over many lattice sites around the impurity. Thus, we expect 
that they will be weakly affected by local modifications of the 
potential near the position of the vacancy. Finally, we have 
not discussed the implications of our results for the magnetic 
properties of graphene planes. An enhancement of the density 
of states over distances which are large compared to the lattice 
spacing, implies that vacancies may lead to the formation of 
extended magnetic moments, enhancing the tendency of the 
system towards ferro or antiferromagnetism. 
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FIG. 3: (color online) Inverse participation ratio for a concentration 
of vacancies of 0.5 % in systems with 10 4 sites. The dashed curve is 
the total DOS for the same dilution. Top: t' = 0; center: t' = 0.1 t; 
bottom: t' = 0.2 1. Only the low energy region is shown. 

the 2D Dirac equation, we find quasi-localized states at the 
Fermi level, if the clean electronic structure shows electron- 
hole symmetry. In the absence of electron-hole symmetry, 
these states are shifted and broadened. In a system with a 
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